#####################################################################
# Normalization of Censorship: Evidence from China
# Replication Files (Experimental Studies, Online Appendices)
#
# Author: Tony Zirui Yang
# Date: January 24, 2024
#
# The lines below reproduce the main Figures and results
# related to the experimental studies reported in the appendices
#####################################################################

rm(list=ls())

# Packages ----

library("stargazer")
library("tidyverse")
library("rstatix")
library("ggpubr")
library("ggplot2")
library("gridExtra")
library("data.table")

# Load the data
Experiment <- fread("Experiment.csv")
# Study 1 data
Study1 <- Experiment[which(Experiment$Study == 1),]
# Study 2 data
Study2 <- Experiment[which(Experiment$Study == 2),]


# --------------------------------------------------
# Table A1
# --------------------------------------------------
# Experiment 1
# Region
a1.1 <- paste(round(table(Study1$Region)/nrow(Study1) * 100, 1), "%", sep = "")
# Gender
a1.2 <- paste(round(rev(table(Study1$Female)/nrow(Study1)) * 100, 1), "%", sep = "")
# Education
# Combining 4 year college and above
Study1$Education4 <- ifelse(Study1$Education == 5, 4, Study1$Education)
a1.3 <- paste(round(table(Study1$Education4)/nrow(Study1) * 100, 1), "%", sep = "")
# Age Group
# Combining age groups into 5 categories
Study1$Age_Group5 <- ifelse(Study1$Age_Group %in% c(2,3), 2, 
                            ifelse(Study1$Age_Group %in% c(4,5), 3,
                                   ifelse(Study1$Age_Group %in% c(6,7), 4, 
                                          ifelse(Study1$Age_Group %in% c(8), 5, 1))))
a1.4 <- paste(round(table(Study1$Age_Group5)/nrow(Study1) * 100, 1), "%", sep = "")
# Income
# Combining 8,000-10,000 and >10,000
Study1$Income4 <- ifelse(Study1$Income == 5, 4, Study1$Income)
a1.5 <- paste(round(table(Study1$Income4)/nrow(Study1) * 100, 1), "%", sep = "")
# Occupation
# Combining unemployed, self-employed, and retired
Study1$Occupation <- ifelse(Study1$Occupation %in% c(11,12,13), 11, Study1$Occupation)
a1.6 <- paste(round(table(Study1$Occupation)/nrow(Study1) * 100, 1), "%", sep = "")
# Location
a1.7 <- paste(round(rev(table(Study1$Urban)/nrow(Study1)) * 100, 1), "%", sep = "")

a1.s1 <- c(a1.1,a1.2,a1.3,a1.4,a1.5,a1.6,a1.7)

# Experiment 2
# Region
a1.8 <- paste(round(table(Study2$Region)/nrow(Study2) * 100, 1), "%", sep = "")
# Gender
a1.9 <- paste(round(rev(table(Study2$Female)/nrow(Study2)) * 100, 1), "%", sep = "")
# Education
a1.10 <- paste(round(table(Study2$Education)/nrow(Study2) * 100, 1), "%", sep = "")
# Age Group
# Combining age groups into 5 categories
Study2$Age_Group5 <- ifelse(Study2$Age_Group %in% c(2,3), 2, 
                            ifelse(Study2$Age_Group %in% c(4,5), 3,
                                   ifelse(Study2$Age_Group %in% c(6,7), 4, 
                                          ifelse(Study2$Age_Group %in% c(8), 5, 1))))
a1.11 <- paste(round(table(Study2$Age_Group5)/nrow(Study2) * 100, 1), "%", sep = "")
# Income
# Combining 8,000-10,000 and >10,000
Study2$Income4 <- ifelse(Study2$Income == 5, 4, Study2$Income)
a1.12 <- paste(round(table(Study2$Income4)/nrow(Study2) * 100, 1), "%", sep = "")
# Occupation
# Combining unemployed and retired
a1.13 <- rep(NA,11)
# Location
a1.14 <- rep(NA,2)

a1.s2 <- c(a1.8,a1.9,a1.10,a1.11,a1.12,a1.13,a1.14)

# China Internet Census ()
# Region
a1.15 <- paste(c(0.462,0.221,0.233,0.084) * 100, "%", sep = "")
# Gender
a1.16 <- paste(c(0.481,0.519) * 100, "%", sep = "")
# Education
a1.17 <- paste(c(0.561,0.238,0.105,0.097) * 100, "%", sep = "")
# Age Group
a1.18 <- paste(c(0.232,0.215,0.208,0.176,0.169) * 100, "%", sep = "")
# Income
a1.19 <- paste(c(0.510,0.215,0.143,0.133) * 100, "%", sep = "")
# Occupation
a1.20 <- paste(c(0.269,0.224,0.080,0.029,0.028,0.060,0.026,0.044,0.042,0.063,0.135) * 100, "%", sep = "")
# Location
a1.21 <- paste(c(0.718,0.282) * 100, "%", sep = "")
a1.c <- c(a1.15,a1.16,a1.17,a1.18,a1.19,a1.20,a1.21)


table.A1 <- cbind(a1.s1,a1.s2,a1.c)
colnames(table.A1) <- c("Study 1", "Study 2", "China Internet Census")
rownames(table.A1) <- c("Region: East", "Region: Central", "Region: West", "Region: Northeast",
                        "Gender: Female", "Gender: Male",
                        "Education: <= Junior High", "Education: Senior High", "Education: 3-year College", "Education: >= 4-year College",
                        "Age: <=19", "Age: 20-29", "Age: 30-39", "Age: 40-49", "Age: >=50",
                        "Income: <=3000", "Income: 3000-5000", "Income: 5000-8000", "Income: >=8000",
                        "Ocuupation: Student", "Ocuupation: Self-employed", "Ocuupation: Corporate employee", "Ocuupation: Corporate management",
                        "Ocuupation: Government employee", "Ocuupation: Professional", "Ocuupation: Manufacturing", "Ocuupation: Service worker",
                        "Ocuupation: Migrant worker", "Ocuupation: Farmer", "Ocuupation: Unemployed & Retired",
                        "Location: Urban", "Location: Rural")
table.A1


# --------------------------------------------------
# Table A2
# --------------------------------------------------
# Covariates for balance check
covariates <- c("Female", "Age_Group", "Education", "Income",
                "Party_Member", "Ideology", "Political_Interest", "Social_Media")
# Study 1
t.tests.1 <- lapply(covariates, function(x) {
  t.test(as.formula(paste(x, "~ Treatment")), data = Study1,
         alternative = c("two.sided", "less", "greater"),
         mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
})
# Study 2
t.tests.2 <- lapply(covariates, function(x) {
  t.test(as.formula(paste(x, "~ Treatment")), data = Study2,
         alternative = c("two.sided", "less", "greater"),
         mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
})
# Combined
t.tests.c <- lapply(covariates, function(x) {
  t.test(as.formula(paste(x, "~ Treatment")), data = Experiment,
         alternative = c("two.sided", "less", "greater"),
         mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
})
# Gender
a2.1 <- round(c(t.tests.1[[1]]$estimate, t.tests.1[[1]]$p.value, 
                t.tests.2[[1]]$estimate, t.tests.2[[1]]$p.value, 
                t.tests.c[[1]]$estimate, t.tests.c[[1]]$p.value),3)
# Age Group
a2.2 <- round(c(t.tests.1[[2]]$estimate, t.tests.1[[2]]$p.value, 
                t.tests.2[[2]]$estimate, t.tests.2[[2]]$p.value, 
                t.tests.c[[2]]$estimate, t.tests.c[[2]]$p.value),3)
# Education
a2.3 <- round(c(t.tests.1[[3]]$estimate, t.tests.1[[3]]$p.value, 
                t.tests.2[[3]]$estimate, t.tests.2[[3]]$p.value, 
                t.tests.c[[3]]$estimate, t.tests.c[[3]]$p.value),3)
# Income
a2.4 <- round(c(t.tests.1[[4]]$estimate, t.tests.1[[4]]$p.value, 
                t.tests.2[[4]]$estimate, t.tests.2[[4]]$p.value, 
                t.tests.c[[4]]$estimate, t.tests.c[[4]]$p.value),3)
# Party Member
a2.5 <- round(c(t.tests.1[[5]]$estimate, t.tests.1[[5]]$p.value, 
                t.tests.2[[5]]$estimate, t.tests.2[[5]]$p.value, 
                t.tests.c[[5]]$estimate, t.tests.c[[5]]$p.value),3)
# Ideology
a2.6 <- round(c(t.tests.1[[6]]$estimate, t.tests.1[[6]]$p.value, 
                t.tests.2[[6]]$estimate, t.tests.2[[6]]$p.value, 
                t.tests.c[[6]]$estimate, t.tests.c[[6]]$p.value),3)
# Political Interest
a2.7 <- round(c(t.tests.1[[7]]$estimate, t.tests.1[[7]]$p.value, 
                t.tests.2[[7]]$estimate, t.tests.2[[7]]$p.value, 
                t.tests.c[[7]]$estimate, t.tests.c[[7]]$p.value),3)
# Social Media Usage
a2.8 <- round(c(t.tests.1[[8]]$estimate, t.tests.1[[8]]$p.value, 
                t.tests.2[[8]]$estimate, t.tests.2[[8]]$p.value, 
                t.tests.c[[8]]$estimate, t.tests.c[[8]]$p.value),3)

table.A2 <- rbind(a2.1,a2.2,a2.3,a2.4,a2.5,a2.6,a2.7,a2.8)
colnames(table.A2) <- c("Study 1: Control", "Study 1: Treated", "Study 1: p",
                        "Study 2: Control", "Study 2: Treated", "Study 2: p",
                        "Combined: Control", "Combined: Treated", "Combined: p")
rownames(table.A2) <- c("Female", "Age Group", "Education", "Income", "Party Member",
                        "Ideology", "Pol Interest", "Social Media")
table.A2


# --------------------------------------------------
# Table A3
# --------------------------------------------------
lm_a3.1 <- lm(Treatment ~ Female + Education + Age_Group + Income + Ideology + Party_Member + 
               Political_Interest + Social_Media, data = Study1)
lm_a3.2 <- lm(Treatment ~ Female + Education + Age_Group + Income + Ideology + Party_Member + 
               Political_Interest + Social_Media, data = Study2)
lm_a3.c <- lm(Treatment ~ Female + Education + Age_Group + Income + Ideology + Party_Member + 
               Political_Interest + Social_Media, data = Experiment)

stargazer(lm_a3.1, lm_a3.2, lm_a3.c, 
          style = "APSR",
          covariate.labels = c("Female", "Education", "Age Group", "Income", "Ideology",
                               "Party Member", "Political Interest", "Social Media Usage"))


# --------------------------------------------------
# Table B1
# --------------------------------------------------
lm_b1.1 <- lm(Censor_Support ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study1)
lm_b1.2 <- lm(Censor_Support ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study2)
lm_b1.3 <- lm(Censor_Support ~ Treatment + Ideology, data = Experiment)
lm_b1.4 <- lm(Censor_Support ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Experiment)

stargazer(lm_b1.1, lm_b1.2, lm_b1.3, lm_b1.4, 
          style = "APSR",
          title = "Support for Censorship Apparatus",
          covariate.labels = c("Treatment", "Female", "Education", "Age Group", "Income", 
                               "Ideology", "Party Member", "Political Interest", "Social Media Usage"))


# --------------------------------------------------
# Table B2
# --------------------------------------------------
lm_b2.1 <- lm(Regime_Overall ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study1)
lm_b2.2 <- lm(Regime_Overall ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study2)
lm_b2.3 <- lm(Regime_Overall ~ Treatment + Ideology, data = Experiment)
lm_b2.4 <- lm(Regime_Overall ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Experiment)

stargazer(lm_b2.1, lm_b2.2, lm_b2.3, lm_b2.4, 
          style = "APSR",
          title = "Regime Support: Overall Satisfaction",
          covariate.labels = c("Treatment", "Female", "Education", "Age Group", "Income", 
                               "Ideology", "Party Member", "Political Interest", "Social Media Usage"))


# --------------------------------------------------
# Table B3
# --------------------------------------------------
lm_b3.1 <- lm(Regime_Central ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study1)
lm_b3.2 <- lm(Regime_Central ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study2)
lm_b3.3 <- lm(Regime_Central ~ Treatment + Ideology, data = Experiment)
lm_b3.4 <- lm(Regime_Central ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Experiment)

stargazer(lm_b3.1, lm_b3.2, lm_b3.3, lm_b3.4, 
          style = "APSR",
          title = "Regime Support: Central Government",
          covariate.labels = c("Treatment", "Female", "Education", "Age Group", "Income", 
                               "Ideology", "Party Member", "Political Interest", "Social Media Usage"))

# --------------------------------------------------
# Table B4
# --------------------------------------------------
lm_b4.1 <- lm(Regime_Local ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study1)
lm_b4.2 <- lm(Regime_Local ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study2)
lm_b4.3 <- lm(Regime_Local ~ Treatment + Ideology, data = Experiment)
lm_b4.4 <- lm(Regime_Local ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Experiment)

stargazer(lm_b4.1, lm_b4.2, lm_b4.3, lm_b4.4, 
          style = "APSR",
          title = "Regime Support: Local Government",
          covariate.labels = c("Treatment", "Female", "Education", "Age Group", "Income", 
                               "Ideology", "Party Member", "Political Interest", "Social Media Usage"))

# --------------------------------------------------
# Table B5
# --------------------------------------------------
lm_b5.1 <- lm(Regime_Protest ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study1)
lm_b5.2 <- lm(Regime_Protest ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Study2)
lm_b5.3 <- lm(Regime_Protest ~ Treatment + Ideology, data = Experiment)
lm_b5.4 <- lm(Regime_Protest ~ Treatment + Female + Education + Age_Group + Income + Ideology + Party_Member + 
                Political_Interest + Social_Media, data = Experiment)

stargazer(lm_b5.1, lm_b5.2, lm_b5.3, lm_b5.4, 
          style = "APSR",
          title = "Willingness to Protest",
          covariate.labels = c("Treatment", "Female", "Education", "Age Group", "Income", 
                               "Ideology", "Party Member", "Political Interest", "Social Media Usage"))


# --------------------------------------------------
# Figure B1
# --------------------------------------------------

# Heterogeneous Treatment Effect on Support for Censorship
Figure.B1.1.1 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.1.2 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Female == 1),])

Figure.B1.1.3 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.1.4 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.1.5 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])

Figure.B1.1.6 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.1.7 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.1.8 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])

Figure.B1.1.9 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.1.10 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.1.11 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])

Figure.B1.1.12 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.1.13 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.1.14 <- lm(Censor_Support ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])


Figure.B1.Terms <- c("Male", "Female", "High school", "3-year college", "4-year college",
           "Age < 30", "Age in 30-40", "Age > 40", "Income < 5K", "Income in 3K-5K", "Income > 8K",
           "Conservative", "Moderate", "Liberal")
Figure.B1.Estimates.1 <- c(coef(summary(Figure.B1.1.1))["Treatment","Estimate"], coef(summary(Figure.B1.1.2))["Treatment","Estimate"],
               coef(summary(Figure.B1.1.3))["Treatment","Estimate"], coef(summary(Figure.B1.1.4))["Treatment","Estimate"],
               coef(summary(Figure.B1.1.5))["Treatment","Estimate"], coef(summary(Figure.B1.1.6))["Treatment","Estimate"],
               coef(summary(Figure.B1.1.7))["Treatment","Estimate"], coef(summary(Figure.B1.1.8))["Treatment","Estimate"],
               coef(summary(Figure.B1.1.9))["Treatment","Estimate"], coef(summary(Figure.B1.1.10))["Treatment","Estimate"],
               coef(summary(Figure.B1.1.11))["Treatment","Estimate"], coef(summary(Figure.B1.1.12))["Treatment","Estimate"],
               coef(summary(Figure.B1.1.13))["Treatment","Estimate"], coef(summary(Figure.B1.1.14))["Treatment","Estimate"])
Figure.B1.Std_Error.1 <- c(coef(summary(Figure.B1.1.1))["Treatment","Std. Error"], coef(summary(Figure.B1.1.2))["Treatment","Std. Error"],
               coef(summary(Figure.B1.1.3))["Treatment","Std. Error"], coef(summary(Figure.B1.1.4))["Treatment","Std. Error"],
               coef(summary(Figure.B1.1.5))["Treatment","Std. Error"], coef(summary(Figure.B1.1.6))["Treatment","Std. Error"],
               coef(summary(Figure.B1.1.7))["Treatment","Std. Error"], coef(summary(Figure.B1.1.8))["Treatment","Std. Error"],
               coef(summary(Figure.B1.1.9))["Treatment","Std. Error"], coef(summary(Figure.B1.1.10))["Treatment","Std. Error"],
               coef(summary(Figure.B1.1.11))["Treatment","Std. Error"], coef(summary(Figure.B1.1.12))["Treatment","Std. Error"],
               coef(summary(Figure.B1.1.13))["Treatment","Std. Error"], coef(summary(Figure.B1.1.14))["Treatment","Std. Error"])
Figure.B1.1 <- data.frame(Figure.B1.Terms, Figure.B1.Estimates.1, Figure.B1.Std_Error.1)

multiplier <- 1.96
for (i in 1:nrow(Figure.B1.1)){
  Figure.B1.1$ymin[i] = Figure.B1.1$Figure.B1.Estimates.1[i] - (multiplier * Figure.B1.1$Figure.B1.Std_Error.1[i])
  Figure.B1.1$ymax[i] = Figure.B1.1$Figure.B1.Estimates.1[i] + (multiplier * Figure.B1.1$Figure.B1.Std_Error.1[i])
}
Figure.B1.1$Figure.B1.Terms <- factor(Figure.B1.1$Figure.B1.Terms, levels=rev(Figure.B1.1$Figure.B1.Terms), ordered=TRUE)

Figure.B1.p1 <-  ggplot(Figure.B1.1, aes(x=Figure.B1.Terms, y=Figure.B1.Estimates.1)) + 
  geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
  geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
  labs(y="Coefficient of Estimates", x="", title="Support for Censorship Apparatus") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw() + 
  theme(axis.text.x = element_text(size = 12, face = "bold"),
        axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme

# Heterogeneous Treatment Effect on Regime Support: Overall Satisfaction
Figure.B1.2.1 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.2.2 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Female == 1),])

Figure.B1.2.3 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.2.4 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.2.5 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])

Figure.B1.2.6 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.2.7 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.2.8 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])

Figure.B1.2.9 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.2.10 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.2.11 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])

Figure.B1.2.12 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.2.13 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.2.14 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])


Figure.B1.2.Estimates <- c(coef(summary(Figure.B1.2.1))["Treatment","Estimate"], coef(summary(Figure.B1.2.2))["Treatment","Estimate"],
               coef(summary(Figure.B1.2.3))["Treatment","Estimate"], coef(summary(Figure.B1.2.4))["Treatment","Estimate"],
               coef(summary(Figure.B1.2.5))["Treatment","Estimate"], coef(summary(Figure.B1.2.6))["Treatment","Estimate"],
               coef(summary(Figure.B1.2.7))["Treatment","Estimate"], coef(summary(Figure.B1.2.8))["Treatment","Estimate"],
               coef(summary(Figure.B1.2.9))["Treatment","Estimate"], coef(summary(Figure.B1.2.10))["Treatment","Estimate"],
               coef(summary(Figure.B1.2.11))["Treatment","Estimate"], coef(summary(Figure.B1.2.12))["Treatment","Estimate"],
               coef(summary(Figure.B1.2.13))["Treatment","Estimate"], coef(summary(Figure.B1.2.14))["Treatment","Estimate"])
Figure.B1.2.Std_Error <- c(coef(summary(Figure.B1.2.1))["Treatment","Std. Error"], coef(summary(Figure.B1.2.2))["Treatment","Std. Error"],
               coef(summary(Figure.B1.2.3))["Treatment","Std. Error"], coef(summary(Figure.B1.2.4))["Treatment","Std. Error"],
               coef(summary(Figure.B1.2.5))["Treatment","Std. Error"], coef(summary(Figure.B1.2.6))["Treatment","Std. Error"],
               coef(summary(Figure.B1.2.7))["Treatment","Std. Error"], coef(summary(Figure.B1.2.8))["Treatment","Std. Error"],
               coef(summary(Figure.B1.2.9))["Treatment","Std. Error"], coef(summary(Figure.B1.2.10))["Treatment","Std. Error"],
               coef(summary(Figure.B1.2.11))["Treatment","Std. Error"], coef(summary(Figure.B1.2.12))["Treatment","Std. Error"],
               coef(summary(Figure.B1.2.13))["Treatment","Std. Error"], coef(summary(Figure.B1.2.14))["Treatment","Std. Error"])
Figure.B1.2 <- data.frame(Figure.B1.Terms, Figure.B1.2.Estimates, Figure.B1.2.Std_Error)

for (i in 1:nrow(Figure.B1.2)){
  Figure.B1.2$ymin[i] = Figure.B1.2$Figure.B1.2.Estimates[i] - (multiplier * Figure.B1.2$Figure.B1.2.Std_Error[i])
  Figure.B1.2$ymax[i] = Figure.B1.2$Figure.B1.2.Estimates[i] + (multiplier * Figure.B1.2$Figure.B1.2.Std_Error[i])
}
Figure.B1.2$Figure.B1.Terms <- factor(Figure.B1.2$Figure.B1.Terms, levels=rev(Figure.B1.2$Figure.B1.Terms), ordered=TRUE)

Figure.B1.p2 <- ggplot(Figure.B1.2, aes(x=Figure.B1.Terms, y=Figure.B1.2.Estimates)) + 
  geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
  geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
  labs(y="Coefficient of Estimates", x="", title="Regime Support: Overall") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw() + 
  theme(axis.text.x = element_text(size = 12, face = "bold"),
        axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme


# Heterogeneous Treatment Effect on Regime Support: Central Government
Figure.B1.3.1 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.3.2 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Female == 1),])

Figure.B1.3.3 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.3.4 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.3.5 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])

Figure.B1.3.6 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.3.7 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.3.8 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])

Figure.B1.3.9 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.3.10 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.3.11 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])

Figure.B1.3.12 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.3.13 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.3.14 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])


Figure.B1.3.Estimates <- c(coef(summary(Figure.B1.3.1))["Treatment","Estimate"], coef(summary(Figure.B1.3.2))["Treatment","Estimate"],
               coef(summary(Figure.B1.3.3))["Treatment","Estimate"], coef(summary(Figure.B1.3.4))["Treatment","Estimate"],
               coef(summary(Figure.B1.3.5))["Treatment","Estimate"], coef(summary(Figure.B1.3.6))["Treatment","Estimate"],
               coef(summary(Figure.B1.3.7))["Treatment","Estimate"], coef(summary(Figure.B1.3.8))["Treatment","Estimate"],
               coef(summary(Figure.B1.3.9))["Treatment","Estimate"], coef(summary(Figure.B1.3.10))["Treatment","Estimate"],
               coef(summary(Figure.B1.3.11))["Treatment","Estimate"], coef(summary(Figure.B1.3.12))["Treatment","Estimate"],
               coef(summary(Figure.B1.3.13))["Treatment","Estimate"], coef(summary(Figure.B1.3.14))["Treatment","Estimate"])
Figure.B1.3.Std_Error <- c(coef(summary(Figure.B1.3.1))["Treatment","Std. Error"], coef(summary(Figure.B1.3.2))["Treatment","Std. Error"],
               coef(summary(Figure.B1.3.3))["Treatment","Std. Error"], coef(summary(Figure.B1.3.4))["Treatment","Std. Error"],
               coef(summary(Figure.B1.3.5))["Treatment","Std. Error"], coef(summary(Figure.B1.3.6))["Treatment","Std. Error"],
               coef(summary(Figure.B1.3.7))["Treatment","Std. Error"], coef(summary(Figure.B1.3.8))["Treatment","Std. Error"],
               coef(summary(Figure.B1.3.9))["Treatment","Std. Error"], coef(summary(Figure.B1.3.10))["Treatment","Std. Error"],
               coef(summary(Figure.B1.3.11))["Treatment","Std. Error"], coef(summary(Figure.B1.3.12))["Treatment","Std. Error"],
               coef(summary(Figure.B1.3.13))["Treatment","Std. Error"], coef(summary(Figure.B1.3.14))["Treatment","Std. Error"])
Figure.B1.3 <- data.frame(Figure.B1.Terms, Figure.B1.3.Estimates, Figure.B1.3.Std_Error)

for (i in 1:nrow(Figure.B1.3)){
  Figure.B1.3$ymin[i] = Figure.B1.3$Figure.B1.3.Estimates[i] - (multiplier * Figure.B1.3$Figure.B1.3.Std_Error[i])
  Figure.B1.3$ymax[i] = Figure.B1.3$Figure.B1.3.Estimates[i] + (multiplier * Figure.B1.3$Figure.B1.3.Std_Error[i])
}
Figure.B1.3$Figure.B1.Terms <- factor(Figure.B1.3$Figure.B1.Terms, levels=rev(Figure.B1.3$Figure.B1.Terms), ordered=TRUE)

Figure.B1.p3 <- ggplot(Figure.B1.3, aes(x=Figure.B1.Terms, y=Figure.B1.3.Estimates)) + 
  geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
  geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
  labs(y="Coefficient of Estimates", x="", title="Regime Support: Central") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw() + 
  theme(axis.text.x = element_text(size = 12, face = "bold"),
        axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme


# Heterogeneous Treatment Effect on Regime Support: Local Government
Figure.B1.4.1 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.4.2 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Female == 1),])

Figure.B1.4.3 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.4.4 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.4.5 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])

Figure.B1.4.6 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.4.7 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.4.8 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])

Figure.B1.4.9 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.4.10 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.4.11 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])

Figure.B1.4.12 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.4.13 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.4.14 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])


Figure.B1.4.Estimates <- c(coef(summary(Figure.B1.4.1))["Treatment","Estimate"], coef(summary(Figure.B1.4.2))["Treatment","Estimate"],
               coef(summary(Figure.B1.4.3))["Treatment","Estimate"], coef(summary(Figure.B1.4.4))["Treatment","Estimate"],
               coef(summary(Figure.B1.4.5))["Treatment","Estimate"], coef(summary(Figure.B1.4.6))["Treatment","Estimate"],
               coef(summary(Figure.B1.4.7))["Treatment","Estimate"], coef(summary(Figure.B1.4.8))["Treatment","Estimate"],
               coef(summary(Figure.B1.4.9))["Treatment","Estimate"], coef(summary(Figure.B1.4.10))["Treatment","Estimate"],
               coef(summary(Figure.B1.4.11))["Treatment","Estimate"], coef(summary(Figure.B1.4.12))["Treatment","Estimate"],
               coef(summary(Figure.B1.4.13))["Treatment","Estimate"], coef(summary(Figure.B1.4.14))["Treatment","Estimate"])
Figure.B1.4.Std_Error <- c(coef(summary(Figure.B1.4.1))["Treatment","Std. Error"], coef(summary(Figure.B1.4.2))["Treatment","Std. Error"],
               coef(summary(Figure.B1.4.3))["Treatment","Std. Error"], coef(summary(Figure.B1.4.4))["Treatment","Std. Error"],
               coef(summary(Figure.B1.4.5))["Treatment","Std. Error"], coef(summary(Figure.B1.4.6))["Treatment","Std. Error"],
               coef(summary(Figure.B1.4.7))["Treatment","Std. Error"], coef(summary(Figure.B1.4.8))["Treatment","Std. Error"],
               coef(summary(Figure.B1.4.9))["Treatment","Std. Error"], coef(summary(Figure.B1.4.10))["Treatment","Std. Error"],
               coef(summary(Figure.B1.4.11))["Treatment","Std. Error"], coef(summary(Figure.B1.4.12))["Treatment","Std. Error"],
               coef(summary(Figure.B1.4.13))["Treatment","Std. Error"], coef(summary(Figure.B1.4.14))["Treatment","Std. Error"])
Figure.B1.4 <- data.frame(Figure.B1.Terms, Figure.B1.4.Estimates, Figure.B1.4.Std_Error)

for (i in 1:nrow(Figure.B1.4)){
  Figure.B1.4$ymin[i] = Figure.B1.4$Figure.B1.4.Estimates[i] - (multiplier * Figure.B1.4$Figure.B1.4.Std_Error[i])
  Figure.B1.4$ymax[i] = Figure.B1.4$Figure.B1.4.Estimates[i] + (multiplier * Figure.B1.4$Figure.B1.4.Std_Error[i])
}
Figure.B1.4$Figure.B1.Terms <- factor(Figure.B1.4$Figure.B1.Terms, levels=rev(Figure.B1.4$Figure.B1.Terms), ordered=TRUE)

Figure.B1.p4 <- ggplot(Figure.B1.4, aes(x=Figure.B1.Terms, y=Figure.B1.4.Estimates)) + 
  geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
  geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
  labs(y="Coefficient of Estimates", x="", title="Regime Support: Local") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw() + 
  theme(axis.text.x = element_text(size = 12, face = "bold"),
        axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme


# Heterogeneous Treatment Effect on Willingness to Protest
Figure.B1.5.1 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.5.2 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Female == 1),])

Figure.B1.5.3 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.5.4 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.5.5 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])

Figure.B1.5.6 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.5.7 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.5.8 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])

Figure.B1.5.9 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.5.10 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.5.11 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])

Figure.B1.5.12 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.5.13 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.5.14 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])


Figure.B1.5.Estimates <- c(coef(summary(Figure.B1.5.1))["Treatment","Estimate"], coef(summary(Figure.B1.5.2))["Treatment","Estimate"],
               coef(summary(Figure.B1.5.3))["Treatment","Estimate"], coef(summary(Figure.B1.5.4))["Treatment","Estimate"],
               coef(summary(Figure.B1.5.5))["Treatment","Estimate"], coef(summary(Figure.B1.5.6))["Treatment","Estimate"],
               coef(summary(Figure.B1.5.7))["Treatment","Estimate"], coef(summary(Figure.B1.5.8))["Treatment","Estimate"],
               coef(summary(Figure.B1.5.9))["Treatment","Estimate"], coef(summary(Figure.B1.5.10))["Treatment","Estimate"],
               coef(summary(Figure.B1.5.11))["Treatment","Estimate"], coef(summary(Figure.B1.5.12))["Treatment","Estimate"],
               coef(summary(Figure.B1.5.13))["Treatment","Estimate"], coef(summary(Figure.B1.5.14))["Treatment","Estimate"])
Figure.B1.5.Std_Error <- c(coef(summary(Figure.B1.5.1))["Treatment","Std. Error"], coef(summary(Figure.B1.5.2))["Treatment","Std. Error"],
               coef(summary(Figure.B1.5.3))["Treatment","Std. Error"], coef(summary(Figure.B1.5.4))["Treatment","Std. Error"],
               coef(summary(Figure.B1.5.5))["Treatment","Std. Error"], coef(summary(Figure.B1.5.6))["Treatment","Std. Error"],
               coef(summary(Figure.B1.5.7))["Treatment","Std. Error"], coef(summary(Figure.B1.5.8))["Treatment","Std. Error"],
               coef(summary(Figure.B1.5.9))["Treatment","Std. Error"], coef(summary(Figure.B1.5.10))["Treatment","Std. Error"],
               coef(summary(Figure.B1.5.11))["Treatment","Std. Error"], coef(summary(Figure.B1.5.12))["Treatment","Std. Error"],
               coef(summary(Figure.B1.5.13))["Treatment","Std. Error"], coef(summary(Figure.B1.5.14))["Treatment","Std. Error"])
Figure.B1.5 <- data.frame(Figure.B1.Terms, Figure.B1.5.Estimates, Figure.B1.5.Std_Error)

for (i in 1:nrow(Figure.B1.5)){
  Figure.B1.5$ymin[i] = Figure.B1.5$Figure.B1.5.Estimates[i] - (multiplier * Figure.B1.5$Figure.B1.5.Std_Error[i])
  Figure.B1.5$ymax[i] = Figure.B1.5$Figure.B1.5.Estimates[i] + (multiplier * Figure.B1.5$Figure.B1.5.Std_Error[i])
}
Figure.B1.5$Figure.B1.Terms <- factor(Figure.B1.5$Figure.B1.Terms, levels=rev(Figure.B1.5$Figure.B1.Terms), ordered=TRUE)

Figure.B1.p5 <- ggplot(Figure.B1.5, aes(x=Figure.B1.Terms, y=Figure.B1.5.Estimates)) + 
  geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
  geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
  labs(y="Coefficient of Estimates", x="", title="Willingness to Protest") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw() + 
  theme(axis.text.x = element_text(size = 12, face = "bold"),
        axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme


Figure.B1 <- grid.arrange(Figure.B1.p1,
                          Figure.B1.p2,
                          Figure.B1.p3,
                          Figure.B1.p4,
                          Figure.B1.p5, nrow = 3)
# Save the figure
ggsave(Figure.B1,
       file = "Figures_Appendix/Figure_B1.pdf",   # The directory you want to save the file in
       width = 9, # The width of the plot in inches
       height = 10) # The height of the plot in inches




# --------------------------------------------------
# Table B6
# --------------------------------------------------
# Outcome variables
formula <- c("Censor_Support ~ ", "Regime_Overall ~ ", "Regime_Central ~ ",
             "Regime_Local ~ ", "Regime_Protest ~ ")
# Study 1 Results (Difference in Means)
lm_list_s1 <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Study1))
names(lm_list_s1) <- c("lm_censor_s1", "lm_overall_s1", "lm_central_s1", "lm_local_s1", "lm_protest_s1")
# Study 2 Results (Difference in Means)
lm_list_s2 <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Study2))
names(lm_list_s2) <- c("lm_censor_s2", "lm_overall_s2", "lm_central_s2", "lm_local_s2", "lm_protest_s2")
# Two Studies Combined Results (Difference in Means)
lm_list_c <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Experiment))
names(lm_list_c) <- c("lm_censor_c", "lm_overall_c", "lm_central_c", "lm_local_c", "lm_protest_c")
# Combining Models into a list
models <- c(lm_list_s1, lm_list_s2, lm_list_c)
# Point Estimates
Estimates <- round(sapply(models, function(x) coef(summary(x))["Treatment","Estimate"]),3)
# p-values
p.value <- paste("[", round(sapply(models, function(x) coef(summary(x))["Treatment","Pr(>|t|)"]),5), "]")
# Adjusting p-values using Benjamini-Hochberg p-value correction
Adj.p.value <- paste("[", round(p.adjust(sapply(models, function(x) coef(summary(x))["Treatment","Pr(>|t|)"]), method = "BH"),5), "]")

table.B6 <- rbind(Estimates[c(1,6,11)], p.value[c(1,6,11)], Adj.p.value[c(1,6,11)],
                  Estimates[c(2,7,12)], p.value[c(2,7,12)], Adj.p.value[c(2,7,12)],
                  Estimates[c(3,8,13)], p.value[c(3,8,13)], Adj.p.value[c(3,8,13)],
                  Estimates[c(4,9,14)], p.value[c(4,9,14)], Adj.p.value[c(4,9,14)],
                  Estimates[c(5,10,15)], p.value[c(5,10,15)], Adj.p.value[c(5,10,15)])
colnames(table.B6) <- c("Study 1", "Study 2", "Combined")
rownames(table.B6) <- c("Support for Censorship Apparatus", "   p-value", "   adjusted p-value",
                        "Overall Satisfaction of China", "   p-value", "   adjusted p-value",
                        "Assessment of Central Governmen", "   p-value", "   adjusted p-value",
                        "Assessment of Local Government", "   p-value", "   adjusted p-value",
                        "Willingness to Protest", "   p-value", "   adjusted p-value")
table.B6


# --------------------------------------------------
# Figure B2
# --------------------------------------------------
# Create a new variable for Explicit Support for Censorship Apparatus (Binary)
Study2 <- Study2 %>% 
  mutate(across(
    c(Censor_Support),
    list(Binary = ~ case_when(
      . %in% c(4,5) ~ 1,
      . %in% c(1,2,3) ~ 0,
      TRUE ~ NA_real_
    ))
  ))
## Group Means for Explicit Censorship Support
Figure.B2 <-  Study2 %>%
  drop_na(Treatment) %>%
  group_by(Treatment) %>%
  get_summary_stats(Censor_Support_Binary, type = "mean_sd") %>%
  mutate(variable = case_when(
    variable == "Censor_Support_Binary" ~ "Explicit Support"
  ),
  Treatment = case_when(
    Treatment == 1 ~ "Treatment",
    Treatment == 0 ~ "Control"
  ),
  # Standard errors for group mean of Explicit Censorship Support
  se = sd / sqrt(n),
  # Confidence interval for group mean Explicit Censorship Support
  ymin = mean - 1.96 * se,
  ymax = mean + 1.96 * se) 

# Calculating Implicit Censorship Support
# Implicit support in the treatment group
lmt <- lm(Censor_Support_Implicit ~ Censor_Support_LE_T, data = Study2[Study2$Treatment == 1,])
# Implicit support in the control group
lmc <- lm(Censor_Support_Implicit ~ Censor_Support_LE_T, data = Study2[Study2$Treatment == 0,])

# Create a data frame for Figure B2
# Data frame for implicit support
Treatment <- c("Control", "Treatment")
Variable <- c(rep("Implicit Support",2))
Estimates <- c(coef(summary(lmc))["Censor_Support_LE_T","Estimate"],
               coef(summary(lmt))["Censor_Support_LE_T","Estimate"])
Std_Error <- c(coef(summary(lmc))["Censor_Support_LE_T","Std. Error"],
               coef(summary(lmt))["Censor_Support_LE_T","Std. Error"])
lower <- c(confint(lmc, level = 0.95)["Censor_Support_LE_T",][1],
           confint(lmt, level = 0.95)["Censor_Support_LE_T",][1])
upper <- c(confint(lmc, level = 0.95)["Censor_Support_LE_T",][2],
           confint(lmt, level = 0.95)["Censor_Support_LE_T",][2])
Figure.B2.2  <- data.frame(Treatment, Variable, Estimates, Std_Error, lower, upper)
# Combining the data frame for explicit support and the data frame for implicit support
Figure.B2  <- Figure.B2[,c(1,2,4,6,7,8)]
colnames(Figure.B2) <- colnames(Figure.B2.2)
Figure.B2 <- rbind(Figure.B2, Figure.B2.2)

# Plot out implicit and explicit support
pd = position_dodge2(width = 0.3, reverse = TRUE)
ggplot(Figure.B2, aes(x=Treatment, y = Estimates, shape = Variable)) + 
  # Point estimates
  geom_point(position = pd, size = 3) + 
  # Confidence intervals
  geom_linerange(aes(x = Treatment, 
                     ymin = lower,
                     ymax = upper),
                 size = 0.5,
                 position = pd, show.legend = FALSE) +
  ggtitle("Support for Censorship Apparatus") +
  labs(y="Proportion of Respondents", x="") +  # Labels
  theme_bw() +
  scale_y_continuous(limits = c(0.2,0.8)) +
  theme(axis.text = element_text(size = 12, face = "bold"),
        axis.title = element_text(size=12),
        axis.text.y = element_text(size = 12),
        legend.title = element_blank()) # Legend label for shape
# Save the figure
ggsave(file = "Figures_Appendix/Figure_B2.pdf",   # The directory you want to save the file in
       width = 5.6, # The width of the plot in inches
       height = 4) # The height of the plot in inches

